All-dielectric polarization-sensitive metasurface for terahertz polarimetric imaging

Terahertz polarimetric imaging, capable of capturing not only intensity profiles but also the polarization states of the incident pattern, is an essential technique with promising applications such as security scans and medical diagnoses. Recently, a novel approach for terahertz imaging has been proposed using a metasurface absorber that converts terahertz light into a temperature profile. However, polarization remains indistinguishable in the imaging process due to the isotropic geometry of the metasurface. To address this issue, this study introduces an all-dielectric, polarization-sensitive metasurface absorber and showcases its suitability for terahertz polarimetric imaging. Optical and thermal simulations confirm that the polarization dependence of our metasurface is translated into the thermal domain, allowing us to distinguish both intensity and polarization states in the incoming image. Additionally, we demonstrate that polarimetric imaging under general, elliptical polarization is attainable. This metasurface facilitates terahertz polarimetric imaging, eliminating the need for complex setups or bulky components, thereby reducing the form factor and enabling widespread use.

Period p = 251.1 μm and height h = 100 μm are fixed for all resonators.

Metasurface for large-scale simulations
Fig. S1.A schematic diagram of the metasurface used for large-scale simulations.

Thermal analysis
The temperature distribution of each individual unit cell in an 8 by 8 arrays, as well as the change in temperature over time, are determined using the stationary equation (Eq.S1) and the time-dependent equation (Eq.S2), respectively.The equations are given as follows: ∇ •  = , (S1) where  is the heat flux, which is expressed in the product of the thermal conductivity of substance and its temperature gradient, i.e.  = −∇, Q is the heat source obtained from the power dissipation density derived from the optical simulations,  is the density of material, and   is the heat capacity.The thermal conductivity, density, and heat capacity of the silicon (Si) resonators, which depend on temperature, are calculated by utilizing the built-in functions of COMSOL Multiphysics.Additionally, the thermal conductivity, density, and heat capacity of polydimethylsiloxane (PDMS) are 0.16 W/m • K, 970 kg/m 3 , and 1460 J/kg • K, respectively, which are the standard material properties of PDMS.
All surfaces of the structure, which consists of Si resonators and the PDMS film, are exposed to the external air, and a Gaussian beam is located in the center of the structure.Each resonator conducts the heat generated by the Gaussian beam to the PDMS underneath, and the resonator and the PDMS eventually dissipate the heat to the environment via thermal radiation and natural convection.The thermal radiation emitted from the surfaces of the structure to the surrounding air is denoted as follows: Note that   is the thermal radiation heat flux,  is the emissivity,  is the Stefan-Boltzmann constant, and Tamb is the temperature of the external air.The boundaries for natural convection are the upper surfaces of each resonator and the upper/lower surfaces of the PDMS film.Because the peripheral side walls of the PDMS film, which is 8-μm -thick, exhibit a negligible airflow compared to those on the upper or lower surfaces of the PMDS film [1], it is justified not to take the contribution of natural convection on the side walls of the PDMS film into account.Natural convection on the side walls of the 100-μm -high resonators is also neglected because the buoyancy of the side wall is significantly weakened, leading to making conduction predominate over the other heat transfer modes [1][2][3].
Heat dissipation via natural convection is governed by Newton's law of cooling, which is expressed as follows: where   is the convective heat flux, ℎ ̅ is the average natural convection coefficient,  ext is the external temperature, and T is the temperature at the surface of the material.The average natural convection coefficients are obtained by two empirical correlations, which are for the cases that natural convection between a hot plate (here, the top surfaces of each resonator and the PDMS) and an upper atmosphere (Eq.S5) and that between a hot plate (here, the bottom surface of the PDMS) and a lower atmosphere (Eq.S6): where  air is the thermal conductivity of air, L is the characteristic length calculated as the area of the PDMS divided by the perimeter of the PDMS, and   is the Rayleigh number.
The average natural convection coefficients, calculated based on Eq.S5, are 4.73 W/ 2 • K for the minimum temperature difference between the external air and the surface of the structure and 11.64 W/ 2 • K for the largest temperature difference.The average natural convection coefficients at the bottom of the PDMS film are calculated based on Eq.S6 and are half of those obtained by Eq.S5.Because the natural convection coefficients of gases generally range from 2-25 W/ 2 • K, the predicted average natural convection coefficients seem to be reasonable [4].where ℎ is the Planck constant,   is the Boltzmann constant,  is the speed of light, and  is the emissivity.We assume that the temperature in the surrounding air is at room temperature (298 K).Then the power captured by an infrared camera, with its focal point positioned at the bottom of the substrate, is computed by using a convolution theorem.The axial full-width halfmaximum of the point spread function is set as 10 μm.By integrating this spectral irradiance with three dimensions-namely, two spatial dimensions and the wavelength-across the wavelength domain, a 2D heat distribution is derived.Finally, the 2D temperature profile is calculated using the equation  IR (, ) = ((, )/) 1 4 , where  is included for solid angle integration.

Robustness of the metasurfaces under fabrication errors
The robustness of the metasurfaces against fabrication errors is presented in this section.
Simulations are conducted for three resonator configurations: x/y-aligned resonator (#1), 45°/-45°-aligned resonator (#2), and RCP/LCP resonator (Fig. S2).For the x/y-aligned resonator, total absorption (A) and absorption ratio between the two resonators (Q0/Q90) are obtained under horizontally polarized light at 1 THz by varying L1 and w1 (Fig. S2A and B).The formula L =  1 +  1 , w =  1 +  1 is used to introduce errors ranging from -10 to 10 μm.These results show that the polarization-sensitive absorption does not degrade significantly upon the fabrication errors.Another noticeable feature is that A and Q0/Q90 exhibit the opposite behavior, that is, the parameters resulting in high A tends to have lower Q0/Q90.Our design is optimized to provide high absorption and strong polarization dependence targeted at 1 THz between this trade-off.Indeed, the changes of A and Q0/Q90 shown in Fig. S2A and B are associated with the frequency shift.As these results at 1 THz does not show spectral response, the spectra of A and Q0/Q90 of the metasurface when dL1 = 10 μm and dw1 = -2 μm are presented in Fig. S2C and   D. This indicates that the absorption peak is blueshifted.While our original metasurface is focused at 1 THz, it can be tuned for other target frequency by adjusting the geometric parameters within the same structure.Similarly, we conduct similar analysis for the 45°/-45°-aligned resonator by changing (L2, w2), (L2, l2), and (w2, l2) (see Fig. S3) and for the RCP/LCP resonator by changing (a, b) and θ (see Fig. S4).Overall, the results show the similar tendency with high sensitivity upon the errors in width (w), which is natural considering that the absolute length of the width is the smallest and therefore the error of -10 to 10 μm corresponds to the large relative errors.Furthermore,

Steady-state temperature rise under various input powers
In the main manuscript, the input power is 0.1 mW for all thermal simulations.For completeness, we simulate the steady-state temperature of the unit cell consisting of 45° and -45°-aligned under various input power (Fig. S5).When applying 45°-polarized incidence, the two resonators show a distinct temperature rise, indicating that our metasurface enables the polarimetric imaging under different light intensities.The thermal response time is conventionally defined as the thermal time constant τ th , indicating that the sensor output signal reaches 63.2% of its the steady-state value [5].we obtain τ th of each resonator by fitting the temperature profiles shown in Fig. 3  where   and   are initial and final temperatures, respectively (Fig. S6).Each resonator approaches the stationary state in different rate and thus has distinct τ th .We define the thermal response time of each resonator as τ th .The thermal response time of each resonator is presented in Table 2 Therefore, the response times of unit cell #1, #2, and #3 are 0.83, 0.73, and 0.85 seconds, respectively.
The three-dimensional (3D) volumetric temperature distribution is acquired through thermal simulation.To transform this 3D distribution into a two-dimensional (2D) infrared image profile, we first calculate the 4D spectral irradiance based on the 3D temperature data.(, , , )

Fig
Fig. S4C and D prove that the fabrication error in angle also does not affect the performance significantly.

Fig. S4 .
Fig. S4.Optical response of the RCP / LCP resonator changing parameters (A, B) a, b and (C, D) θ.For clarification, unit cells at five different points are shown in (E).

Table 2 .
. Thermal response time of each resonator.